StatQuest生物统计学专题 - PCA

PCA的简明含义PCA的内在算法PCA在R中的实现

本周开始，将会学习机器学习的内容，此后会接触很多相关的概念：分类、聚类、降维，有监督、无监督等等。

本周主题是PCA(principal component analysis,主成分分析)，PCA是经典的降维方法之一，在多维变量场景中使用的频次非常高。

先来看一下PCA的含义。

PCA的简明含义

直白来讲，PCA就是把需要通过成百上千个变量描述的数据归结为只需要少数几个变量（经典的是2个）来描述数据。

也就是“两句话”概括“千言万语”。

以学生及学生的各科成绩为例来说明：

假定每个学生共有6科成绩，那么为了描述每一个学生成绩就需要6个变量。但是如果进行PCA分析之后，就只用2个变量（PC1，PC2）就可以描述每一个学生的成绩。

# 原始数据
student        subject1    subject2    subject3    subject4    subject5    subject6
Jim            90          80          70          60          50          40  
Joe            75          90          60          30          40          60
Lily        80          95          45          88          89          30

# PCA后
student        PC1     PC2
Jim            -1.0    1.6
Joe            -1.2    -1.5
Lily        2.2     -0.1

如下图所示，原本需要用6个变量描述的学生成绩只需要2个变量就可以很好的描述。

“把六句话概括为两句话”。

StatQuest-XIII-1

此时两个主成分可以解释100%的原数据波动，此时的碎石图（scree plot）如下：

从左自右分别是PC1、PC2、PC3，三者所能解释的数据波动分别为62%、38%、0%。

StatQuest-XIII-2

上述过程的R操作如下：

#构造数据
student <- c(paste0(rep('subject',6),seq(1,6)))
student
Jim <- c(90,80,70,60,50,40)
Joe <- c(75,90,60,30,40,60)
Lily <- c(80,95,45,88,89,30)
stu.data <- t(data.frame(Jim=Jim,Joe=Joe,Lily=Lily))
colnames(stu.data) <- student
#PCA分析
pca.data <- prcomp(stu.data,scale. = T)
#PCA图
plot(pca.data$x[,1],pca.data$x[,2],xlim = c(-2.0,3.0),xlab = 'PC1',ylab = 'PC2')
text(pca.data$x[,1]+0.3,pca.data$x[,2],rownames(pca.data$x))
#PCA的各成分R2
SS <- pca.data$sdev^2
pca.R2 <- round(SS/sum(SS)*100,1)
barplot(pca.R2,xlim=c(0,6),ylim=c(0,80))
text(seq(1,3),pca.R2+5,paste0(pca.R2,"%"))

本例中，原数据只有6个变量，因此并没有特别显示出PCA的降维功效。

假如是小鼠的单细胞RNA-seq数据，每个细胞有10,000个基因Read Count值，那么使用PCA将原本的10,000个变量降维为2个变量来描述就会很有意义。

一般而言，PCA分析之后会给出一个PCA图，而这个图往往会显示出某种分群的特性，就如小鼠的单细胞RNA-seq来说，小鼠的不同细胞（如神经细胞、上皮细胞）大概率在PCA图上分别聚团。

如下图所示，不同类别的细胞分别聚团。

StatQuest-XIII-3

需要明确的是PCA分析后得出的主成分（PC1、PC2等等）都是没有现实的意义的，你很难说他们代表什么。

不过他们都是降维变量（上例中的6个科目成绩）的线性组合。比如，

PC2和PC3同理。

# 各个主成分的载入因子 
          PC1   PC2   PC3
subject1 -0.07  0.66 -0.11
subject2  0.37 -0.46  0.51
subject3 -0.46  0.30  0.53
subject4  0.46  0.31 -0.24
subject5  0.51  0.09 -0.16
subject6 -0.41 -0.40 -0.60

# R实现
round(pca.data$rotation,2)

PC1的载入因子绝对值最大是subject5，因此可以粗略的说PC1是subject5的代表。同理PC2受到subject1和subject2影响较大，可以认为PC2是subject1和subject2的综合代表。

可以通过这种形式对每个主成分有一个粗略的概念，然而大概率这些主成分都没有现实意义，它们只是一种数学衡量。

如果想将数据降维为一些有现实意义的因子，那么应该考虑使用探索性因子分析（Exploratory Factor Analysis，EFA）。

PCA的内在算法

以小鼠的RNA-seq数据为例，简单起见，先看2个gene的情形：

StatQuest-XIII-4

注意：本例是对基因Read值进行PCA降维分析。不过上表中基因是行变量，这种变量安排和R中是相反的，在R中，列向量是进行PCA降维分析的变量，行变量是PCA图上的各个点。

PCA的算法如下：

1. 寻找中心点，并平移至中心点；

2. 计算PC1；

   PC1是一条过原点的最佳拟合曲线，所有数据点与其的距离平方和最小。

3. 计算PC2；

   PC2与PC1相垂直。

4. 以PC1和PC2为标准坐标轴，绘制PCA图。

下面分别拆解上述步骤：

1 寻找并平移中心点

如下图所示，先对Gene1求平均值，在对Gene2求平均值，于是就可以获得所有数据的中心点，然后将坐标系原点平移至此中心点。

StatQuest-XIII-5

2 计算PC1

计算PC1就是为了获得一条过原点的直线，所有数据点距离此直线的距离平方和最小。以一个点为例，也就是下图中1图的距离a最小，由于每个数据点距离原点的距离是固定的，所以使得距离a最小，也就是使得距离b最大。

也就是最大。

最终可以求得最佳拟合直线如下图3所示，假定此直线斜率为0.25，也就是说"PC1由4份Gene1和1份Gene2构成"，计算方法是“奇异值分解”。

由于奇异值分解时，斜边c是标准化为1的，所以标准化之后，PC1由”0.97份Gene1和0.242份Gene2构成"

StatQuest-XIII-6

3 计算PC2

在PCA中，各个主成分之间是相互垂直的，在这里PC2和PC1是垂直的，也就是PC2的斜率为4，经过标准化之后，"PC2由-0.242份Gene1和0.97份Gene2构成"。

于是PC1和PC2都已经找到，分别是下图中的红色虚线和蓝色虚线坐标轴。

StatQuest-XIII-7

4 绘制PCA图

于是所有的数据点都可以转换为坐标（PC1，PC2），以PC1和PC2为坐标轴即可绘制出相应的PCA图：

StatQuest-XIII-8

2个Gene的PCA过程如上，3个Gene的PCA同理，先获得中心点，然后对所有点进行拟合获得PC1，然后在垂直于PC2的线中，找到最佳拟合线，最后PC3是垂直于PC1和PC2的。

拟合的要求就是过中心点，同时点与拟合直线的距离平方和最小()，也就是使最大。

有几个样本就会有几个主成分，且主成分的重要性依次下降，PC1>PC2>PC3。

PCA在R中的实现

1 首先构造虚拟数据

构造的数据的行是小鼠，列是基因，也就是说列是用于求解主成分的变量。

data.matrix <- matrix(nrow=100, ncol=10)
colnames(data.matrix) <- c(paste0("wt", 1:5),  paste0("ko", 1:5))
rownames(data.matrix) <- paste0("gene", 1:100)
for (i in 1:100) {
  wt.values <- rpois(5, lambda=sample(x=10:1000, size=1))
  ko.values <- rpois(5, lambda=sample(x=10:1000, size=1))

  data.matrix[i,] <- c(wt.values, ko.values)
}
data.matrix <- t(data.matrix)

2 函数prcomp用于实施PCA分析

pca <- prcomp(data.matrix, scale=TRUE)

3 绘制PCA图

pca中有一个变量x，它的行是小鼠，列是不同的主成分。

只绘制PC1和PC2的PCA图：

StatQuest-XIII-9

4 绘制碎石图

pca中有一个变量sdev，它是各个主成分的标准差。

pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)

barplot(pca.var.per, ylim=c(0,105),main="Scree Plot",
        xlab="Principal Component", ylab="Percent Variation")
text(1:length(pca.var.per)+seq(-0.2,1.5,length.out = length(pca.var.per)),
     pca.var.per+5,
     paste0(pca.var.per,'%'))

如下图所示，

StatQuest-XIII-10

5 获得对PC1贡献最大的前10个基因

PC1是由100个Gene的线性方程组成的，各个Gene的系数在pca的rotataion变量中。

loading_scores <- pca$rotation[,1]
gene_scores <- abs(loading_scores) 
gene_score_ranked <- sort(gene_scores, decreasing=TRUE)
top_10_genes <- names(gene_score_ranked[1:10])

top_10_genes 
#Result："gene69" "gene81" "gene8"  "gene21" "gene97" "gene38" "gene88" "gene64" "gene33" "gene96"

本周的PCA分享就到这里，PCA其实是通过“奇异值分解”来完成的，下周就来深入探讨一下PCA的“奇异值分解”的R实现过程。

专题以往文章

参考资料

StatQuest课程：https://statquest.org/video-index/

---

 猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。